In [17]:
set.seed(999)
options(scipen = 9)
options(warn = -1)
source("./environment/libraries.R")
knitr::opts_chunk$set(fig.height = 12, fig.width = 18, fig.dpi = 300)
knitr::opts_chunk$set(warning = FALSE)
All packages loaded successfully
In [18]:
name <- "Kenya_E1"
dataset <- read.csv(paste0("./test/", name, "_processed.csv"))
dataset_c <- data.frame(dataset[7:ncol(dataset)], row.names = dataset$Serial)
proj_coord <- data.frame("Easting" = dataset$Easting, "Northing" = dataset$Northing)
xy_coord <- data.frame("Longitude" = dataset$Longitude, "Latitude" = dataset$Latitude)
dataset_c_closed <- cbind(dataset_c, "Res" = 100 - rowSums(dataset_c))
head(dataset_c_closed)
dataset_spdf <- cbind(proj_coord, dataset_c_closed)
rownames(dataset_spdf) <- rownames(dataset_c)
structures_geojson_path <- file.path("./data", paste0(name, "_topo_lines.geojson"))
structures <- st_read(structures_geojson_path, quiet = TRUE)
structures_utm <- st_transform(structures, crs = 32637) # Transform to UTM Zone 35S (WGS84)
output_dirs <- "./test/ck/"
if (!dir.exists(output_dirs)) {
dir.create(output_dirs, recursive = TRUE)
}
| P | Ca | Ti | Mn | Fe | Ni | Zn | Rb | Sr | Y | Zr | Ba | Mg | Al | Si | K | Res | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | |
| 1 | 0.2731 | 3.4795 | 0.3309 | 0.0603 | 2.7248 | 0.0021 | 0.0033 | 0.0035 | 0.0808 | 0.0017 | 0.012300000 | 0.0933 | 1.227527 | 5.454678 | 21.96007 | 1.971416 | 62.32071 |
| 2 | 0.3835 | 3.9587 | 0.2828 | 0.0566 | 2.4054 | 0.0017 | 0.0031 | 0.0037 | 0.0839 | 0.0016 | 0.004100000 | 0.0955 | 1.234100 | 4.319073 | 18.23019 | 2.227139 | 66.70890 |
| 3 | 0.4158 | 4.0720 | 0.2265 | 0.0617 | 2.1634 | 0.0011 | 0.0033 | 0.0034 | 0.0921 | 0.0012 | 0.004547633 | 0.0755 | 1.316831 | 4.052775 | 17.19447 | 2.489834 | 67.82554 |
| 4 | 0.3704 | 3.9589 | 0.2902 | 0.0601 | 2.6374 | 0.0022 | 0.0035 | 0.0030 | 0.0837 | 0.0021 | 0.010200000 | 0.0880 | 1.384247 | 4.668478 | 18.72138 | 2.425924 | 65.29028 |
| 5 | 0.2981 | 3.3979 | 0.2984 | 0.0707 | 2.7003 | 0.0026 | 0.0030 | 0.0035 | 0.0797 | 0.0017 | 0.015200000 | 0.1005 | 1.280048 | 4.962708 | 20.89937 | 2.486514 | 63.39976 |
| 6 | 0.4373 | 4.0080 | 0.1808 | 0.0537 | 1.9270 | 0.0010 | 0.0029 | 0.0030 | 0.0860 | 0.0012 | 0.003527636 | 0.0583 | 1.304530 | 3.233936 | 13.66489 | 2.973143 | 72.06078 |
In [19]:
source("./utils/functions/create_quick_map.R")
options(repr.plot.width = 10, repr.plot.height = 10)
dbscan_result <- dbscan(dataset_spdf[, c("Easting", "Northing")], eps = 10, minPts = 5)
Area <- factor(dbscan_result$cluster)
dataset$Area <- Area
create_quick_map(dataset, structures, group_data = "Area")
dataset_spdf$Area <- Area
dataset_spdf <- dataset_spdf[dataset_spdf$Area == 2, ]
In [20]:
create_spatial_objects <- function(data_spdf, transformation = "none") {
coords <- data_spdf[, c("Easting", "Northing")]
crs_string <- "+proj=utm +zone=37 +north +datum=WGS84"
if (transformation == "ilr") {
spatial_data <- as.data.frame(ilr(data_spdf[, -c(1, 2, ncol(data_spdf))]))
} else if (transformation == "clr") {
spatial_data <- as.data.frame(clr(data_spdf[, -c(1, 2, ncol(data_spdf))]))
} else {
spatial_data <- data_spdf[, -c(1, 2, ncol(data_spdf))]
}
return(SpatialPointsDataFrame(coords = coords, data = spatial_data,
proj4string = CRS(crs_string)))
}
spdf_comp <- create_spatial_objects(dataset_spdf, "none")
spdf_ilr <- create_spatial_objects(dataset_spdf, "ilr")
spdf_clr <- create_spatial_objects(dataset_spdf, "clr")
source("./utils/functions/lag_distance_from_spdf.R")
source("./utils/functions/site_diagonal_from_spdf.R")
lag_dist <- lag_distance_from_spdf(spdf_ilr)
site_diag <- site_diagonal_from_spdf(spdf_ilr)
In [21]:
source("./utils/functions/geostatistical_analysis.R")
source("./utils/functions/validate_geostatistical_predictions.R")
In [22]:
results_ilr_omni <- geostatistical_analysis(spdf_ilr, "ILR_Omnidirectional", lag_dist, site_diag, directional = FALSE)
Selected models:
psill range kappa n
Nug 0.005636354 0.000000 0 136
Sph 0.013267064 4.879795 0 69
Gau 0.006745542 3.073336 0 38
Exp 0.018271788 1.914671 0 29
Linear Model of Coregionalization found. Good.
[using ordinary cokriging]
Saved 17 maps to: ./test/ck/ilr-omnidirectional/maps/
In [23]:
validate_predictions(spdf_ilr, results_ilr_omni$ck, results_ilr_omni$g)
V1 V2 V3 V4 V5
ME 0.0001022348 -0.0062020027 -0.0013457263 -0.0008333475 0.0007015358
MSE 0.008724575 0.072485737 0.018183152 0.013407582 0.037057207
V6 V7 V8 V9 V10
ME 0.0015601790 0.0031828449 0.0014674274 -0.0015878812 -0.0082124051
MSE 0.008616293 0.021176479 0.027666382 0.018340586 0.103614886
V11 V12 V13 V14 V15
ME 0.0014876697 0.0015760092 0.0003796410 0.0011482589 0.0017138209
MSE 0.023532552 0.012736559 0.021056677 0.027255291 0.011694586
V16
ME 0.0007488469
MSE 0.010896231
- $mv_results
A data.frame: 1 × 3 Accuracy Precision Goodness <dbl> <dbl> <dbl> 0.85 0.64 0.82 - $var_results
A data.frame: 16 × 4 Variable Accuracy Precision Goodness <chr> <dbl> <dbl> <dbl> V1 0.95 0.86 0.93 V2 0.85 0.78 0.89 V3 0.90 0.87 0.93 V4 0.95 0.79 0.89 V5 0.95 0.93 0.97 V6 0.85 0.95 0.97 V7 0.90 0.90 0.95 V8 0.95 0.83 0.92 V9 0.95 0.87 0.93 V10 0.90 0.90 0.95 V11 0.90 0.86 0.93 V12 0.95 0.84 0.92 V13 0.95 0.78 0.89 V14 0.95 0.80 0.90 V15 0.95 0.83 0.92 V16 0.95 0.91 0.95
In [24]:
results_ilr_dir <- geostatistical_analysis(spdf_ilr, "ILR_Directional", lag_dist, site_diag, directional = FALSE)
Selected models:
psill range kappa n
Nug 0.005636354 0.000000 0 136
Sph 0.013267064 4.879795 0 69
Gau 0.006745542 3.073336 0 38
Exp 0.018271788 1.914671 0 29
Linear Model of Coregionalization found. Good.
[using ordinary cokriging]
Saved 17 maps to: ./test/ck/ilr-directional/maps/
In [25]:
validate_predictions(spdf_ilr, results_ilr_dir$ck, results_ilr_dir$g)
V1 V2 V3 V4 V5
ME -0.0005531453 0.0048508733 0.0024681402 0.0026194095 0.0043327897
MSE 0.008576134 0.078193724 0.018882226 0.013745652 0.038442514
V6 V7 V8 V9 V10
ME 0.0016470742 0.0017199244 -0.0029722507 0.0016671107 0.0004749151
MSE 0.008548680 0.022354672 0.028832000 0.017899969 0.101678939
V11 V12 V13 V14 V15
ME 0.0005650498 0.0016746456 0.0042720624 0.0037208622 0.0015264638
MSE 0.024453100 0.013906207 0.019529877 0.026214458 0.012438465
V16
ME -0.0006585043
MSE 0.011510467
- $mv_results
A data.frame: 1 × 3 Accuracy Precision Goodness <dbl> <dbl> <dbl> 0.85 0.64 0.82 - $var_results
A data.frame: 16 × 4 Variable Accuracy Precision Goodness <chr> <dbl> <dbl> <dbl> V1 0.95 0.86 0.93 V2 0.85 0.78 0.89 V3 0.90 0.87 0.93 V4 0.95 0.79 0.89 V5 0.95 0.93 0.97 V6 0.85 0.95 0.97 V7 0.90 0.90 0.95 V8 0.95 0.83 0.92 V9 0.95 0.87 0.93 V10 0.90 0.90 0.95 V11 0.90 0.86 0.93 V12 0.95 0.84 0.92 V13 0.95 0.78 0.89 V14 0.95 0.80 0.90 V15 0.95 0.83 0.92 V16 0.95 0.91 0.95
In [26]:
par(bg = "white")
options(repr.plot.width = 15, repr.plot.height = 12)
# Compute MAF analysis
g_temp <- create_gstat_from_spdf(spdf_ilr, method = "ordinary")
v_omni <- variogram(g_temp, width = lag_dist / 2, cutoff = site_diag / 3, cross = TRUE)
maf_result <- Maf(spdf_ilr@data, vg = v_omni)
source("./utils/functions/quick_pca_screeplot.R")
quick_pca_screeplot(maf_result)
maf_dataset <- maf_result$scores[, 1:9]
spdf_maf <- SpatialPointsDataFrame(coords = dataset_spdf[, c("Easting", "Northing")],
data = as.data.frame(maf_dataset),
proj4string = CRS("+proj=utm +zone=37 +north +datum=WGS84"))
results_maf_omni <- geostatistical_analysis(spdf_maf, "MAF_Omnidirectional", lag_dist, site_diag,
directional = FALSE, maf_result = maf_result,
resolution = 2)
Selected models:
psill range kappa n
Nug 0.4284907 0.000000 0 45
Exp 0.4890494 1.198881 0 24
Sph 0.7692933 3.407950 0 21
Linear Model of Coregionalization found. Good.
[using ordinary cokriging]
Saved 17 maps to: ./test/ck/maf-omnidirectional/maps/
In [27]:
validate_predictions(spdf_maf, results_maf_omni$ck, results_maf_omni$g)
maf1 maf2 maf3 maf4 maf5
ME 0.007983487 0.009709664 0.015301634 0.012281329 0.007186930
MSE 0.9761022 1.0474878 0.9136258 0.9094230 0.8691944
maf6 maf7 maf8 maf9
ME -0.008298787 0.021281565 -0.007840009 -0.004253346
MSE 0.8925875 0.6713064 0.7971482 0.9827595
- $mv_results
A data.frame: 1 × 3 Accuracy Precision Goodness <dbl> <dbl> <dbl> 0.95 0.54 0.77 - $var_results
A data.frame: 9 × 4 Variable Accuracy Precision Goodness <chr> <dbl> <dbl> <dbl> maf1 0.90 0.93 0.97 maf2 0.95 0.79 0.89 maf3 0.95 0.91 0.95 maf4 0.95 0.83 0.91 maf5 0.65 0.94 0.96 maf6 0.95 0.89 0.94 maf7 0.90 0.81 0.91 maf8 0.95 0.58 0.79 maf9 0.90 0.92 0.96
In [28]:
g_temp2 <- create_gstat_from_spdf(spdf_ilr, method = "ordinary")
v_dir <- variogram(g_temp2, width = lag_dist / 2, cutoff = site_diag / 3,
alpha = c(0, 45, 90, 135), tol.hor = 22.5, cross = FALSE)
maf_result2 <- Maf(spdf_ilr@data, vg = v_dir)
quick_pca_screeplot(maf_result2)
maf_dataset2 <- maf_result2$scores[, 1:10]
spdf_maf2 <- SpatialPointsDataFrame(coords = dataset_spdf[, c("Easting", "Northing")],
data = as.data.frame(maf_dataset2),
proj4string = CRS("+proj=utm +zone=37 +north +datum=WGS84"))
results_maf_dir <- geostatistical_analysis(spdf_maf2, "MAF_Directional", lag_dist, site_diag,
directional = TRUE, maf_result = maf_result2,
resolution = 2)
Directional approach not fully implemented yet - using omnidirectional as fallback
Selected models:
psill range kappa n
Nug 0.5618877 0.000000 0 55
Sph 0.4160659 3.029266 0 34
Gau 0.9858814 1.388666 0 13
Exp 0.5556532 1.361615 0 8
Linear Model of Coregionalization found. Good.
[using ordinary cokriging]
Saved 17 maps to: ./test/ck/maf-directional/maps/
In [29]:
validate_predictions(spdf_maf2, results_maf_dir$ck, results_maf_dir$g)
maf1 maf2 maf3 maf4 maf5
ME 0.008505926 0.010039250 0.004414058 -0.007083254 -0.010341529
MSE 0.9789648 0.9862980 1.0416295 0.9492662 0.8339911
maf6 maf7 maf8 maf9 maf10
ME -0.010617249 0.021731741 0.006564275 -0.001078079 0.014045939
MSE 0.9829563 0.9810363 0.6349502 0.9448917 0.7574230
- $mv_results
A data.frame: 1 × 3 Accuracy Precision Goodness <dbl> <dbl> <dbl> 0.95 0.49 0.74 - $var_results
A data.frame: 10 × 4 Variable Accuracy Precision Goodness <chr> <dbl> <dbl> <dbl> maf1 0.95 0.93 0.96 maf2 0.80 0.96 0.98 maf3 0.95 0.83 0.92 maf4 0.95 0.90 0.95 maf5 0.95 0.76 0.88 maf6 0.90 0.89 0.94 maf7 0.95 0.92 0.96 maf8 0.95 0.77 0.88 maf9 0.95 0.65 0.83 maf10 0.95 0.81 0.91
In [30]:
source("./utils/functions/create_ck_maplist_from_list.R")
ck_files <- list.files("./test/ck/", full.names = TRUE, pattern = paste0(name, ".*\\.rds$"), recursive = TRUE)
ck_list <- lapply(ck_files, readRDS)
if (length(ck_list) > 0) {
cat("Available cokriging results:", length(ck_list), "\n")
cat("Files:", basename(ck_files), "\n")
if (length(ck_list) > 1) {
cat("Multiple methods completed. Individual results saved in method subdirectories\n")
}
} else {
cat("No cokriging results found\n")
}
cat("Geostatistical workflow completed successfully!\n")
Available cokriging results: 4 Files: Kenya_E1_CK2.rds Kenya_E1_CK2.rds Kenya_E1_CK2.rds Kenya_E1_CK2.rds Multiple methods completed. Individual results saved in method subdirectories Geostatistical workflow completed successfully!